function dX = Dynamics( t , X )

global m S CD g rou0 Hre RE

  col = 0    ;
  
  v   = X(1) ;
  the = X(2) ;
  x   = X(3) ;
  y   = X(4) ;
  
  rou = rou0 ; % * exp( ( RE - y ) / Hre ) ;
  
  dv  = - rou * v^2 * S * CD / 2 / m - g * sin( the ) ;
  
dthe  = ( col - g * cos( the ) ) / v ;

  dx  =   v * cos( the ) ;
  
  dy  =   v * sin( the ) ;

  dX  = [ dv ; dthe ; dx ; dy ] ;
  
end